interneurons <- readRDS("../../data/interneurons_kb.rds")
DefaultAssay(interneurons) <- "RNA"
interneurons = FindVariableFeatures(interneurons, selection.method = "vst", nfeatures = 10000)
all.genes = rownames(interneurons)
interneurons = ScaleData(interneurons)
interneurons = RunPCA(interneurons, features = VariableFeatures(object = interneurons))
interneurons = RunUMAP(interneurons, reduction = "pca", dims = 1:25)
interneurons= FindNeighbors(interneurons, reduction = "pca", dims = 1:25)
interneurons = FindClusters(interneurons, resolution = .9)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 19414
## Number of edges: 639483
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8362
## Number of communities: 23
## Elapsed time: 3 seconds
DimPlot(interneurons, reduction = "umap", label = TRUE, repel = TRUE,pt.size = 1)
plot_colors = read_delim("../../data/palettes/interneuron_timepoint_pal.txt",
delim = "/n",col_names = F) %>% pull(X1)
time_points <- unique(interneurons@meta.data["timepoint"]) %>% pull(timepoint)
col_n <-set_names(plot_colors,levels(time_points))
#,order = rev(levels(time_points))
DimPlot(interneurons,
reduction = "umap",
repel = TRUE,
group.by = "timepoint",
pt.size = 1,
cols = col_n,order = rev(levels(time_points)))+
ggtitle(label = "")+ theme(legend.position = "none")
pre.critical <- levels(time_points)[1:5]
# Add new timepoint column for pre and post critical period
interneurons@meta.data <- interneurons@meta.data %>%
mutate(timepoint2 = ifelse(timepoint %in% pre.critical,
"E14-P5",
"P7-Adult"))
DimPlot(interneurons, reduction = "umap", label = TRUE, repel = TRUE,pt.size = 1, group.by = "timepoint2")
DimPlot(interneurons, reduction = "umap", label = TRUE, repel = TRUE, split.by = "timepoint2",pt.size = 1)
saveRDS(interneurons,"../../data/interneurons_labeled.RDS")
DefaultAssay(interneurons) <- "RNA"
all.markers <- FindAllMarkers(interneurons,only.pos = TRUE,logfc.threshold = 0.25)
sig.markers <- all.markers %>%
filter(p_val_adj < 0.05)
sig.markers %>%
filter(cluster == 6) %>%
arrange(desc(avg_log2FC)) %>%
head(20) %>%
kbl() %>%
kable_minimal()
|
|
p_val
|
avg_log2FC
|
pct.1
|
pct.2
|
p_val_adj
|
cluster
|
gene
|
|
Adarb2
|
0
|
13.983543
|
0.824
|
0.530
|
0
|
6
|
Adarb2
|
|
mt-Rnr2
|
0
|
11.129866
|
0.970
|
0.854
|
0
|
6
|
mt-Rnr2
|
|
Hnrnpa2b1
|
0
|
8.444997
|
0.942
|
0.807
|
0
|
6
|
Hnrnpa2b1
|
|
Rpl21
|
0
|
7.708343
|
0.856
|
0.392
|
0
|
6
|
Rpl21
|
|
Tpt1
|
0
|
5.285972
|
0.943
|
0.824
|
0
|
6
|
Tpt1
|
|
Filip1l
|
0
|
4.589018
|
0.653
|
0.127
|
0
|
6
|
Filip1l
|
|
Rpl81
|
0
|
4.502508
|
0.957
|
0.842
|
0
|
6
|
Rpl8
|
|
Cpne4
|
0
|
4.200330
|
0.289
|
0.154
|
0
|
6
|
Cpne4
|
|
Atrx
|
0
|
3.748112
|
0.868
|
0.594
|
0
|
6
|
Atrx
|
|
Ccdc88a
|
0
|
3.693222
|
0.886
|
0.649
|
0
|
6
|
Ccdc88a
|
|
Fubp1
|
0
|
3.588860
|
0.797
|
0.655
|
0
|
6
|
Fubp1
|
|
Rps241
|
0
|
3.502384
|
0.965
|
0.854
|
0
|
6
|
Rps24
|
|
Srrm2
|
0
|
3.416876
|
0.936
|
0.716
|
0
|
6
|
Srrm2
|
|
Rpl39
|
0
|
3.393382
|
0.905
|
0.798
|
0
|
6
|
Rpl39
|
|
Thra
|
0
|
3.275159
|
0.780
|
0.543
|
0
|
6
|
Thra
|
|
Rpsa-ps10
|
0
|
3.041653
|
0.371
|
0.083
|
0
|
6
|
Rpsa-ps10
|
|
Dync1i11
|
0
|
3.034032
|
0.603
|
0.230
|
0
|
6
|
Dync1i1
|
|
Dclk11
|
0
|
2.961476
|
0.669
|
0.412
|
0
|
6
|
Dclk1
|
|
Sfpq1
|
0
|
2.950321
|
0.847
|
0.589
|
0
|
6
|
Sfpq
|
|
Nedd4
|
0
|
2.869791
|
0.824
|
0.476
|
0
|
6
|
Nedd4
|
write_csv(sig.markers, file = "sig_markers_1.csv")
# genes in the P7 clusters that have a log2 fold change >= 1
P7_genes = sig.markers %>%
filter(cluster ==6 & avg_log2FC>=1)
# Convert gene symbols to ENSEMBL
geneconvert = bitr(rownames(P7_genes), fromType="SYMBOL", toType="ENSEMBL", OrgDb="org.Mm.eg.db")
# Convert to MGI IDs
mgi_ids <- select(org.Mm.eg.db, geneconvert$ENSEMBL, "MGI","ENSEMBL") %>% pull(MGI)
# Write MGI IDs for genewalk analysis
mgi_ids %>% as_tibble() %>% write_delim(file = "P7_genes.txt",col_names = FALSE)
go_enrich <- enrichGO(gene = geneconvert$ENSEMBL,
OrgDb = org.Mm.eg.db,
keyType = "ENSEMBL",
pvalueCutoff = 0.05,
ont = "BP",
pAdjustMethod ="BH",readable = TRUE)
# clusterProfiler::cnetplot(go_enrich,
# showCategory = 15,
# layout = "gem",
# cex_gene=2,
# cex_category=2)
go_plot <- barplot(go_enrich,showCategory = 15,font.size = 24)
#scale_y_discrete(expand=c(0.2, 1))
go_plot$data <- go_plot$data %>% tail(-2)
go_plot
# dotplot(go_enrich,
# split="ONTOLOGY",showCategory=15)+
# facet_grid(ONTOLOGY ~ .,scales = "free")+
# scale_y_discrete(expand=c(0.1, 0))
table(Idents(interneurons)) %>% kbl(col.names = c("Cluster","Cell Count")) %>% kable_minimal()
|
Cluster
|
Cell Count
|
|
0
|
2104
|
|
1
|
1882
|
|
2
|
1754
|
|
3
|
1442
|
|
4
|
1348
|
|
5
|
1285
|
|
6
|
1206
|
|
7
|
1149
|
|
8
|
1009
|
|
9
|
894
|
|
10
|
857
|
|
11
|
838
|
|
12
|
723
|
|
13
|
659
|
|
14
|
572
|
|
15
|
492
|
|
16
|
401
|
|
17
|
273
|
|
18
|
256
|
|
19
|
112
|
|
20
|
93
|
|
21
|
35
|
|
22
|
30
|
DimHeatmap(interneurons, dims = 1:15, cells = 500, balanced = TRUE)
DimPlot(interneurons, reduction = "pca", group.by = "timepoint2",pt.size = 1)
DimPlot(interneurons, reduction = "pca", group.by = "timepoint2",pt.size = 1,dims=c(2,3))
DimPlot(interneurons, reduction = "pca", group.by = "timepoint2",pt.size = 1,dims=c(3,4))
DimPlot(interneurons, reduction = "pca", group.by = "timepoint2",pt.size = 1,dims=c(4,5))
DimPlot(interneurons, reduction = "pca", group.by = "timepoint2",pt.size = 1,dims=c(5,6))
DimPlot(interneurons, reduction = "pca", group.by = "timepoint2",pt.size = 1,dims=c(5,7))
DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1,dims=c(7,8))
interneurons[["percent.olf"]] <- PercentageFeatureSet(interneurons, pattern = "^Olfr")
interneurons[["percent.vmn"]] <- PercentageFeatureSet(interneurons, pattern = "^Vmn")
FeaturePlot(interneurons,features = "percent.olf",label = TRUE, repel = TRUE, min.cutoff = "q10", max.cutoff = "q90",pt.size = 1)
FeaturePlot(interneurons,features = "percent.vmn",label = TRUE, repel = TRUE, min.cutoff = "q10", max.cutoff = "q90",pt.size = 1)
